function [bmus,errors] = neuralGasBmusCosine(M, D, selectedBmus)

dataVectors = size(D,1);
units = size(M,1);

% selected bmus
if ((nargin < 3) || (isempty(selectedBmus))), 
  selectedBmus = 1; 
end

bmus = zeros(dataVectors,length(selectedBmus));
errors = bmus;

for i=1:size(D,1)
    x = D(i,:);
    onesM = ones(size(M,1),1);
    xNorm = sqrt(x*x');
    Dist = mod_distance(acos(sum((M.*x(onesM,:)),2)./(sqrt(sum(M.^2,2)).*xNorm)));
  
    if all(selectedBmus==1), [error bmu] = min(Dist); else [error bmu] = sort(Dist); end
    
    if units==1, bmus(i,:) = 1; else bmus(i,:) = bmu(selectedBmus,:)'; end
    errors(i,:) = Q(selectedBmus,:);
end

